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ABSTRACT 
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This  paper  develops  a special  partitioning  method  for  solving  LP 
problems  with  embedded  network  structure.  These  problems  include  many  of 
the  large-scale  LP  problems  of  practical  importance,  particularly  in  the 
fields  of  energy,  scheduling,  and  distribution.  The  special  partitioning 
method,  called  the  simplex  special  ordered  network  (SON)  procedure, 
applies  to  LP  problems  that  contain  both  non-network  rows  and  non-network 
columns,  with  no  restriction  on  the  form  of  the  rows  and  columns  that  do 
not  exhibit  a network  structure.  These  LP/embedded  network  problems 
include  as  a special  case  multi-commodity  network  problems  and  constrained 
network  problems  previously  treated  in  the  literature,  by  simultaneously 
encompassing  both  side  constraints  and  side  activities. 

The  simplex  SON  procedure  solves  these  problems  by  exploiting  the 
network  topology  of  the  basis,  whose  properties  are  characterized  by  means 
of  a specially  constructed  master  basis  tree.  A set  of  fundamental  ex- 
change rules  are  developed  which  Identify  admissible  ways  to  modify  the 
master  basis  tree,  and  hence  the  composition  of  the  partitioned  basis  in- 
verse. The  simplex  SON  method  implements  these  rules  by  a set  of  stream- 
lined labeling  algorithms,  while  maintaining  the  network  portion  of  the  basis 
as  large  as  possible,  thereby  accelerating  the  basis  exchange  step.  As  a 
result,  LP/embedded  network  problems  can  be  solved  with  less  computer 
storage  and  significantly  greater  efficiency  than  available  from  standard 
linear  programming  methods. 

Preliminary  computational  results  are  reported  for  an  all-FORTRAN 


I 


i 


implementation  o£  the  simplex/ SON  algorithm  called  PNET/LP.  The  test 
problems  are  real-world  models  of  physical  distribution  and  scheduling 
systems.  PNET/LP  has  solved  problems  with  6200  rows  and  22,000  columns 
in  less  than  6 minutes,  counting  all  1/0,  on  an  AMDAHL  V-6  with  a FORTRAN 
G compiler. 


1.0  INTRODUCTION 


The  dramatic  successes  of  the  past  several  years  In  solving  pure 
network  problems  [2,  3,  5,  6,  8,  19,  30,  41]  have  motivated  consideration 
of  methods  for  solving  more  general  linear  programming  (LP)  problems 
with  embedded  network  structure.  For  example.  In  the  realm  of  pure  net- 
works (capacitated  minimum  cost  flow  problems) , the  computational  study 
[21]  demonstrates  that  special  purpose  network  computer  codes  are  150-200 
times  faster  than  the  state-of-the-art  LP  code,  APEX  III.  Subsequent 
studies  of  "singly  constrained"  networks  (LP  problems  consisting  of  a net- 
work plus  one  additional  side  constraint)  demonstrated  that  specialized 
methods  also  yield  substantial  computational  advantages  for  problems  that 
do  not  exhibit  pure  network  structures,  but  which  are  "almost  networks." 

The  papers  [18,  20,  31]  show  that  these  problems  can  be  solved  25-50  times 
faster  than  APEX-III.  Many  practical  LP  problems,  however,  contain  em- 
bedded networks  with  multiple  side  constraints  and  multiple  side  variables, 
and  so  it  is  extremely  Important  to  determine  whether  an  efficient  special- 
ized method  can  be  developed  for  these  problems.  The  purpose  of  this 
paper  is  to  describe  such  a method,  called  the  simplex  special  ordered 
network  (SON)  algorithm.  The  simplex  SON  algorithm  is  a primal  basis 
partitioning  method  that  employs  special  updating  and  labeling  procedures 
to  accelerate  computations  involving  the  network-LP  Interface.  A pre- 
liminary FORTRAN  Implementation  of  this  method  solves  real-world  physical 
distribution  models  25-50  times  faster  than  APEX-III,  confirming  that  it  is 
possible  to  create  a marriage  of  network  and  LP  methodology  that  has  ad- 
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vantages  for  more  general  problems  with  embedded  network  structure. 

For  definitional  purposes,  we  refer  to  an  LP/embedded  network  problem 

as  a capacitated  or  uncapacitated  linear  program  in  which  the  coefficient 

matrix  A can  be  characterized  as  follows.  The  A matrix  contains  m + q rows 

and  n + p columns  which  are  ordered  and  scaled  such  that  each  column  of 

the  m x n submatrix  A , consisting  of  the  first  m rows  and  n columns  of 

mn 

A,  has  at  most  one  +1,  one  -1,  and  zeros  elsewhere.  A major  portion  of 

the  LP  literature  has  been  devoted  to  problems  in  which:  (a)  m * n =*  0 

(standard  LP  problems);  (b)  p ■ 0 (multicommodity  networks  and  constrained 

network  problems);  (c)  p - 0 and  the  submatrix  A contains  only  one  non- 

ton 

zero  entry  per  column  (LP/generalized  upper  bounding  (GUB)  problems);  (d) 
p = q * 0 (pure  network  problems). 

The  success  with  special  classes  of  LP/embedded  network  problems, 
already  noted,  has  led  to  speculations  [10,  18,  21]  that  good  results  can 
also  be  obtained  by  extending  these  ideas  for  problems  where  p and  q are 
less  than  n and  m,  respectively.  Motivation  for  such  an  extension  that 
leads  to  a highly  efficient  implementation  has  come  from  a number  of  the 
major  practitioners  of  linear  programming  including  Harvey  Greenberg, 

Milton  Gutterman,  Alan  Goldman,  and  A.  C.  Williams.  In  addition,  the 
members  of  SHARE,  and  a number  of  industrial  and  governmental  agencies 
that  have  large-scale  LP/embedded  network  problems,  have  strongly  stressed 
the  need  for  such  a method.  Several  of  these  individuals  and  groups  have 
urged  us  to  undertake  such  a development,  leading  to  the  results  reported 
in  this  paper.  One  of  the  important  applications  to  which  the  simplex  SON 
method  is  relevant  is  the  national  energy  model  PIES,  developed  by  the 
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Federal  Energy  Agency  [24].  In  the  PIES  model,  q » 0,  ■ • 2500,  n ■ 4400, 
and  p ■ 4500. 

In  general,  it  Is  our  experience  that  most  large-scale  LP  problems 
involving  production  scheduling,  physical  distribution,  facility  location, 

, personnel  assignment,  or  personnel  promotion  contain  a large  embedded 

network  component,  sometimes  consisting  of  several  smaller  embedded  net- 
works. Coupling  constraints  (q  > 0)  arise,  for  example,  from  economies  of 
scale,  limitation  on  the  total  number  of  promotions,  capacity  restrictions 
on  modes  of  transportation  (e.g.,  pipelines,  barges),  limitations  on 
shared  resources,  multiple  criteria,  or  from  combining  the  outputs  of  sub- 
divisions to  meet  overall  demands.  Coupling  columns  (p  > 0)  arise  from 
activities  which  involve  different  time  periods  (e.g.,  storage),  production 
alternatives  (e.g.,  refinery  activities),  or  which  involve  different  sub- 
divisions (e.g.,  assembly).  For  example,  Agrico  Chemical  Fertilizer 
Company  has  physical  distribution  and  facility  location  problems  where 
p ■ 0,  m * 6200,  n “ 22,000,  and  q * 20. 


1.1  History  of  Methods  for  Solving  Special  Classes  of  LP/Embedded  Network 
Problems 


There  are  two  basic  approaches  which  have  been  employed  to  develop 
specialized  techniques  for  the  above  special  classes  (cases  b,  c,  and  d) 
of  LP/embedded  network  problems — decomposition  and  partitioning  methods. 
Decomposition  approaches  are  further  characterized  as  price-directive 
or  resource-directive.  The  papers  [3,  11,  12,  13,  16,  43,  44,  45,  47] 
give  variations  of  price-directive  decomposition  and  the  papers  [3,  14, 
34,  40,  42]  are  resource-directive  decomposition. 
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Partitioning  approaches  can  be  divided  into  partitioning  block 
diagonal  structured  linear  problems  or  general  partitioning  to  exploit 
embedded  substructure  within  general  linear  problems.  Because  they 


share  the  same  basic  principles,  we  will  briefly  review  the  literature 
of  both. 

The  general  idea  of  partitioning  block  diagonal  structured  programs 
was  originally  proposed  by  Dantzig  [14a]  and  in  a slightly  different 
setting,  by  Charnes  and  Cooper  [9,  vol.  2].  Later  Bennett  [7],  Hartman 
and  Lasdon  [26],  Heesterman  [28],  Kaul  [33],  and  Weber  and  White  [46] 
independently  developed  primal  solution  procedures  of  more  general 
scope.  The  paper  by  Hartman  and  Lasdon  [26]  contains  an  excellent  de- 
scription of  this  approach  and  procedures  for  handling  the  working  in- 
verse. Further,  their  paper  contains  computational  results  of  such  an 
algorithm.  Grigorladis  and  Ritter  proposed  a dual  method  [25]. 

Dantzig  and  Van  Slyke  [15]  then  proposed  their  well  known  "GUB" 
specialization  of  the  primal  simplex  algorithm  for  the  case  where  each 
block  contains  only  one  row.  The  general  block  diagonal  procedures  were 
further  refined  and  specialized  for  multicommodity  network  problems  by 
Saigal  [38].  This  specialization  involves  carrying  a working  basis  in- 
verse whose  size  need  not  exceed  the  number  of  saturated  arcs.  Hartman 
and  Lasdon  [27]  developed  efficient  procedures  for  updating  this  working 
basis  inverse.  However,  neither  Saigal  nor  Hartman  and  Lasdon  discuss 
how  this  procedure  may  be  efficiently  implemented.  Maier  [36]  refined 
their  procedures  and  initiated  implementation  discussions.  Kennlngton, 
et  al  [1,  29,  30]  streamlined  the  implementation  procedures  of  Maier 
and  conducted  extensive  computational  testing. 
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Charnes  and  Cooper  [9,  vol.  2]  originally  proposed  partitioning  to 

exploit  embedded  substructure  within  general  linear  problems  in  their 

Double  Reverse  Method.  Bakes  [4]  independently  proposed  an  analogous 

algorithm  when  p ■ 0.  Klingman  and  Russell  [35]  developed  additional 

specializations  for  exploiting  pure  network  substructure  and  Hultz  and 

Klingman  [31,  32]  for  generalized  network  substructure  (i.e.,  where  each 

column  of  A may  have  at  most  two  arbitrary  non-zero  coefficients), 
mn 

These  papers  initiated  the  first  in-depth  discussion  of  implementation 
techniques.  The  paper  by  Glover,  Karney,  Klingman,  and  Russell  [20] 
presents  the  first  computational  results  of  such  an  algorithm  for  pure 
network  substructure  and  the  paper  by  Hultz  and  Klingman  [31]  presents 
the  first  computational  results  for  generalized  network  substructure. 

The  work  by  McBride  [37]  and  Graves  and  McBride  [23]  on  Factorization 
subsequently  redeveloped  and  refined  the  general  procedures  of  [9,  14a]. 
Further,  McBride's  dissertation  [37]  discussed  specialization  of  these 
procedures  for  exploiting  pure  network  substructure. 

1.2  Form  of  the  Simplex  SON  Method  for  LP/Embedded  Network  Problems 

The  simplex  SON  method  constitutes  a highly  efficient  way  to  modify 

and  implement  the  steps  of  the  primal  simplex  algorithm  for  the  completely 

general  case  of  embedded  pure  network  problems  (where  p > 0 and  q > 0). 

The  efficiency  is  the  direct  result  of  exploiting  the  pure  network  portion 

A of  the  coefficient  matrix  and  the  network-LP  Interface  by  special 
mn 

labeling  and  updating  procedures. 

The  starting  point  for  the  algorithm,  following  the  natural  design  of 
partitioning  methods,  is  to  subdivide  the  coefficient  matrix  into  network 
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and  non-network  components.  By  reference  to  this  subdivision,  a basis  in- 
verse compactification  procedure  is  employed  that  maintains  a working  basis 
inverse , V \ whose  dimension  equals  m + q less  the  rank  of  the  basic  sub- 
columns of  A associated  with  A . The  size  of  V ^ therefore  varies 

mn 

dynamically.  This  is  one  of  the  major  features  of  this  algorithm  that 
distinguishes  it  from  partitioning  algorithms  designed  for  constrained 
networks  [31,  32,  35]. 

The  basic  variables  not  associated  with  V ^ are  stored  in  a special 
graph  form  called  the  master  basis  tree.  The  development  of  the  master 
basis  tree  and  the  procedures  for  using  it  to  efficiently  replace  arith- 
metic operations  constitute  the  principal  contributions  of  this  paper. 

We  show  that  the  operations  normally  performed  by  using  the  full  basis 
inverse  can  instead  be  performed  by  special  labeling  and  graph  traversal 
techniques  [5]  applied  to  the  master  basis  tree  and  its  interface  with  V . 
The  organization  of  the  simplex  SON  method  maintains  the  network  portion 


of  the  basis  as  large  as  possible  at  each  iteration,  thereby  enabling 
these  labeling  and  list  processing  procedures  to  operate  on  a maximally 
dimensioned  part  of  the  basis.  This  in  turn  minimizes  the  size  of  V 


The  resulting  advantages  over  the  standard  LP  implementation  approach 
are  several.  First,  the  graph  traversal  operations  reduce  both  the  amount 
of  work  needed  to  perform  the  algorithmic  steps  and  the  amount  of  computer 
memory  required  to  store  essential  data.  Second,  the  algorithm  orients  the 
execution  of  operations  in  a manner  that  is  best  suited  to  the  design  of 
computers  (making  extensive  use  of  linked  list  structures,  pointers,  and 
logical  operaions  in  place  of  arithmetic  operations.)  Third,  the  method 
is  less  susceptible  to  round-off  error  and  numerical  inaccuracy.  Most  of 
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the  operations  are  performed  using  original  problem  data  and  only  the 
residual  portion  of  the  basis  inverse  associated  with  V ^ is  subject  to 
the  slower  customary  updating,  with  its  greater  attendant  susceptibility 
to  round-off  error  and  numerical  inaccuracy.  (The  graph  traversal  pro- 
cedures also  automatically  eliminate  checking  or  performing  arithmetic 
operations  on  zero  elements.) 

2.0  PROBLEM  NOTATION 

The  LP/embedded  network  problem  and  its  dual  may  be  stated  mathematic- 
ally as  follows : 


Primal 


Minimize 
subject  to: 


c x + c x 
n n p p 


A x + A x 
mn  n mp  p 


A x + A x 
qn  n qp  p 

0 < x — u 
n n 


0 < x < u 
P P 


b 

m 


b 

q 


Dual 


Maximize 
subject  to 


wb  +wb  -yu  - ¥u 
mm  q q n p 


wA  + w A - y < c 
m mn  q qn  n 


wA  +wA  - ¥ £.  c 
m mp  q qp  p 


w , w — unrestricted 
m q 


Y,  ¥ ^ 0 


(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 

(9) 


(10) 
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where  A is  (m  x n),  A is  (in  x p),  A is  (q  x n) , and  A is  (q  x p). 
mn  mp  qn  qp 

The  remaining  vectors  are  conformable  vectors  whose  subscripts  indicate 
their  dimensionality. 

Each  column  of  the  matrix  A contains  at  most  one  +1,  one  -1,  and 

mn 

zeros  elsewhere.  Thus  A corresponds  to  t pure  network  problem.  For 

mn 

this  reason  the  (2)  portion  of  the  LP/embedded  network  problem  will  be 

referred  to  as  node  constraints  or  simply  nodes.  The  x^  variables  will  be 

referred  to  as  arc  variables  or  simply  arcs.  The  arcs  will  be  further 

classified  as  ordinary  arcs,  which  have  exactly  two  non-zero  entries  in 

A , and  as  slack  arcs,  which  have  exactly  one  non-zero  entry  in  A 
mn  1 mn 

In  graph  terminology,  a simple  graph  G(V,E)  is  a finite  set  of  ver- 
tices V and  a finite  set  of  edges  E connecting  the  vertices.  Each  element 
of  E is  identified  with  an  unordered  pair  of  distinct  elements  of  V. 
Schematically,  each  edge  connects  two  distinct  vertices,  which  are  then 
considered  to  be  adjacent.  If  the  edge  set  E is  expanded  to  contain  edges 
which  have  both  endpoints  incident  on  the  same  vertex  or  multiple  edges 
connecting  the  same  two  vertices  (parallel  edges),  then  G(V,E)  is  called  a 
general  graph.  If  each  edge  of  the  general  graph  has  an  implied  direction 
then  the  graph  is  sometimes  called  a digraph  or  directed  graph,  the  vertices 
are  referred  to  as  nodes,  and  the  edges  are  referred  to  as  arcs. 

The  underlying  pure  network  problem  A defines  one  or  more  connected 

mn 

digraphs  as  follows.  Each  row  of  A corresponds  to  a node  and  each  column 

mn 

to  an  arc.  The  -1  entry  in  a column  indicates  the  node  where  the  arc  begins 
(from  node)  and  the  +1  entry  in  a column  Indicates  the  node  where  the  arc 
ends  (to  node).  If  a column  has  only  one  non-zero  entry  (a  slack  arc) 
the  endpoints  of  the  arc  are  incident  on  the  same  node.  This  represents- 
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tion  of  A may  consist  of  several  disjoint  connected  digraphs, 
mn 


Because  each  arc  is  associated  with  a variable,  it  has  lover  and 
upper  bounds  and  an  objective  function  coefficient. 


3.0  BASIS  STRUCTURE 


Using  the  standard  bounded  variable  simplex  algorithm,  a basis  B 
for  the  LP/embedded  network  problem  is  a matrix  composed  of  a linearly 
independent  set  of  column  vectors  selected  from  the  coefficient  matrix 


A - 


p 

1 

• 

1 

A 

1 

1 

A 

mn 

1 

mp 

A 

I 

1 

A 

L qn 

1 

1 

1 

qp  J 

The  variables  associated  with  the  column  vectors  of  B are  considered  to 


be  basic  variables  xg  and  all  others  are  non-basic  variables  x^  at  their 


lower  or  upper  bound. 

Without  loss  of  generality,  it  will  be  assumed  that  A has  full  row 
rank.  Any  basis  B for  the  LP/embedded  network  problem  will,  therefore, 
be  a nonsingular  matrix  of  order  (m  + q)  x (m  + q).  Clearly,  any  basis 
matrix  B may  be  partitioned  as  follows: 


B - 


XB1 

j XB2 

’ B11 

L!±l 

. B21 

j B22  . 

(11) 


where  B, , is  a nonsingular  submatrix  of  A 

11  _ _ mn 


Thus  the  basic  variables 
mn 

fBnl 

x_,  associated  with  the  _ columns  are  exclusively  arc  variables. 

21  r^i 

The  basic  variables  x^  associated  with  the  | | columns  may  also  con- 

tain arc  variables. 


(xB1  and  x^  are  written  above  their  associated 
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components  of  B In  (11).) 

Based  on  the  indicated  partitioning  of  (11),  the  basis  inverse 
B * may  be  written  as  follows: 


,-l 


B"ll  + B'u  B12  V_1  B21  B'll 


- V_1  B21  "'ll 


B1l  B12  V 


-1 


.-1 


(12) 


where  V - - Bn  B u B12. 

The  motivation  for  this  way  of  partitioning  (11)  is  to  factor  out 
the  submatrix  B^  of  and  thereby  exploit  its  inherent  triangularity  by 
viewing  and  storing  B^  as  a digraph.  This  handling  of  B^  has  several  ad- 
vantages: (1)  the  graph  contains  only  the  nonzero  components  of  B^,  (2) 

any  operations  involving  B^  may  be  performed  by  traversing  the  associated 
digraph  using  appropriately  designed  labeling  techniques,  (3)  since  B^  con- 
tains only  original  problem  data,  numerical  errors  are  reduced. 

The  only  matrices  required  to  generate  the  basis  inverse,  as  seen  from 
(12),  are  B,  B ^ and  V The  efficiency  of  generating  the  needed  compo- 
nents of  B ^ in  performing  the  steps  of  the  simplex  method  depends  on  the  size 
and  techniques  used  to  store  B^,  the  labeling  procedures  used  with  B^  to 
eliminate  matrix  multiplications  involving  B J , and  the  procedures  used  to 
maintain  the  position  of  B. 

The  developments  of  the  following  sections  show  how  to  handle  these 
considerations  effectively. 


3.1  Graph  Representation  of  Bn 

As  in  the  case  of  A , B,  defines  a digraph.  The  structure  of  this 
ran  11 

digraph  is  a set  of  disjoint  spanning  trees,  each  of  which  is  augmented 
by  an  arc  whose  from  node  and  to  node  are  the  same  node  in  B Such  an 

lit 
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arc  of  B^  will  be  called  a simple  loop.  (This  structure  implies  that 
is  a block  diagonal  matrix  and  each  block  is  triangular.  The  blocks 
each  consist  of  a spanning  tree  plus  a simple  loop.  Thus,  each  block  is 
a quasi-tree.)  This  structural  property  of  B^  follows  directly  from 
the  fact  that  B^  is  a square  nonsingular  matrix  and  that  B^  is  a sub- 
matrix of  A . Thus  each  column  of  B, , contains  at  most  two  non-zero 
mn  11 

unit  entries  of  the  opposite  sign. 


A column  of  B^  corresponding  to  a simple  loop  may  be  of  two  types 

according  to  whether  the  associated  column  of  A is  a slack  arc  or  an 

ordinary  arc.  If  the  A column  is  an  ordinary  arc,  then  the  partltion- 

mn 

ing  of  B has  split  this  arc  between  B^  and  B2J — that  is,  one  of  its 
unit  entries  lies  in  B^  and  the  other  in  B^^.  Because  of  this,  the  algo- 
rithm stores  and  uses  B^  by  keeping  a larger  digraph  than  the  one  corres- 
ponding to  B^.  This  larger  digraph,  called  the  master  basis  tree,  con- 
tains every  node  in  A plus  another  node  called  the  master  root ; thus, 

mn 

it  always  contains  m + 1 nodes  and  m arcs.  The  nodes  of  this  tree  that 

correspond  to  rows  of  since  they  are  external  to  the  nonsingular 

network  structure  of  B^,  are  called  externalized  roots  (ER' s). 

The  master  basis  tree  contains  all  of  the  ordinary  arcs  in  B^. 

Simple  loops  in  B^  are  contained  in  the  master  tree  in  a modified  form. 

If  the  simple  loop  is  a slack  arc  of  A , then  the  simple  loop  is  re- 

mn 

placed  by  an  arc  between  the  master  root  and  its  unique  node.  If  the 

simple  loop  is  an  ordinary  arc  in  A it  is  replaced  by  an  arc  between 

mn 

its  nodes  in  A . Such  area,  thus,  join  ER's  to  nodes  in  B, To  complete 
mn  li 

the  master  basis  tree  each  ER  is  connected  to  the  master  root  by  an  ex- 
ternalized arc  (EA).  Figure  1 graphically  depicts  a master  basis  tree. 
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It  is  important  to  stress  that  the  master  basis  tree  is  a conceptual 
scheme  designed  to  allow  the  simplex/SON  algorithm  to  efficiently  main- 
tain the  partitioning  of  B while  keeping  the  portion  at  maximum  size 
during  each  iteration.  This  construction  should  not  be  confused  with  the 
simple  model  device  sometimes  employed  in  pure  network  settings,  where  a 
pseudo  root  is  added  for  the  purpose  of  giving  each  slack  arc  two  end- 
points. The  connections  represented  by  the  master  basis  tree  Include 
both  network  and  "extra  network"  structures  (mediated  by  externalized 
roots  and  arcs),  and  the  rules  for  operating  on  these  structures  are  of  a 
very  special  type. 
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To  understand  the  function  of  these  rules,  it  should  be  noted  that 
a multiplicity  of  cases  must  be  considered  when  executing  a basis  ex- 
change step.  For  example,  if  an  arc  variable  is  to  be  added  to  the 
basis,  or  transferred  from  xg2  to  x^,  the  endpoint (s)  of  the  arc  may 
consist  of  (1)  two  nodes  in  a single  block  (quasi-tree)  in  B^,  (2)  nodes 
in  different  blocks  (quasi-trees)  in  B^,  (3)  a node  in  and  a node 
in  B (4)  two  nodes  in  B21,  (5)  a single  node  in  B^,  (6)  a single  node 

in  B„,  (where  cases  (5)  and  (6)  apply  to  slack  arcs  of  A ).  A similar 
21  nm 

multiplicity  of  cases  applies  to  removing  an  arc  from  the  basis,  and  the 
combinations  that  result  from  both  adding  and  removing  arcs  are  still  more 
numerous. 

The  use  of  the  master  basis  tree  permits  all  of  these  cases  to  be 
unified  in  a particularly  convenient  fashion.  The  rules  characterizing 
the  conditions  for  adding  and  deleting  arcs,  and  specifying  the  appro- 
priate restructuring  of  the  master  basis  tree,  are  as  follows. 

3. 2 Fundamental  Exchange  Rules 

1.  An  arc  of  xfl2  can  admissably  be  added  to  the  B^  portion  of  the 
basis,  without  deleting  another,  if  an  only  if  its  loop  in  the  master  basis 
tree  contains  at  least  one  EA.  (Such  a loop  can  contain  at  most  two  EA's.) 
The  updated  form  of  the  master  basis  tree  then  occurs  in  the  following 
manner:  (a)  Add  the  new  arc  and  drop  any  EA  from  the  loop,  (b)  Change 
the  status  of  the  ER  formerly  met  by  the  dropped  EA  to  that  of  an  ordinary 
node,  transferring  its  row  from  the  B^  to  the  B^  portion  of  the  basis. 
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2.  An  arc  can  be  deleted  from  B,  (removing  a component  of  x ,) 

11  cl 

without  adding  another  as  follows:  (a)  Identify  the  node  of  the 
selected  arc  that  is  farthest  from  the  master  root.  (b)  Change  this 
node  into  an  ER  node  by  moving  this  node  to  and  attaching  it  to  the 
master  root  by  a newly  created  EA.  At  the  same  time  delete  the  selected 
arc  from  B^. 

3.  An  arc  can  be  added  to  B^  and  another  simultaneously  removed 
from  B^i  as  follows:  (a)  If  the  loop  in  the  master  basis  tree  created  by 
the  added  arc  includes  the  arc  to  be  dropped,  then  the  exchange  step  is 
handled  exactly  as  an  exchange  step  of  an  ordinary  network  basis.  (Thus 
no  EA's  are  added  or  dropped,  and  no  nodes  alter  their  status  as  ordinary 
nodes  or  ER  nodes.)  (b)  If  the  loop  in  the  master  basis  tree  created  by 
the  added  arc  does  not  include  the  arc  to  be  dropped,  then  the  exchange 
may  be  performed  as  a two-part  process  that  applies  the  preceding  rules 

1 and  2 in  either  order  (as  long  as  the  exchange  is  valid) . 

4.  B^  and  B2^  can  be  restructured,  without  adding  or  deleting  basis 
arcs  xfi^,  by  an  exchange  step  that  drops  any  EA  and  adds  another  EA  to  any 
node  of  the  isolated  tree  (excluding  the  master  root)  created  by  dropping 
the  first.  This  step  is  accomplished  by  interchanging  the  ER  status  and 
ordinary  node  status  of  two  nodes  which  swaps  their  corresponding  rows  in 

*n  snd  B2i- 

It  should  be  remarked  that  the  EA's  have  a special  interpretation  in 

these  rules.  Since  the  master  basis  tree  spans  all  nodes  of  A , and 

mn 

always  contains  the  same  total  number  of  arcs  (Including  EA's),  the  number 
of  EA's  corresponds  to  the  number  of  non-arc  variables  in  the  basis  (ele- 
ments of  Xg2>  that  are  required  to  give  the  basis  full  row  rank  for  the 


k 
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A portion  of  the  problem, 
mn 

It  may  also  be  observed  that  the  A portion  of  each  column  of 

mn 

(the  portion  associated  with  the  ER's)  contains  at  most  one  non- 
zero entry.  In  particular,  this  partial  column  is  the  zero  vector  if 
its  associated  arc  (element  of  x^)  does  not  meet  an  ER,  and  otherwise 
contains  a -1  entry  if  the  from  node  of  the  arc  meets  an  ER  and  a +1 

entry  if  the  to  node  of  the  arc  meets  an  ER.  No  slack  arcs  of  A can 

mn 

contain  entries  in  or  e*8e  £heir  columns  would  be  zero  vectors. 

However,  it  is  possible  for  slack  arcs  of  A to  have  entries  in 

mn 

B_„,  and  also  for  ordinary  arcs  of  A to  have  both  of  their  non-zero 
ii  mn 

entries  in  (meeting  two  ER's  at  their  endpoints).  In  this  case 
such  a slack  arc  or  ordinary  arc,  when  added  to  the  master  basis  tree, 
creates  a loop  that  Includes  an  EA,  and  thus  can  be  moved  from  x^  to 
by  the  Fundamental  Exchange  Rules. 

The  validity  of  the  Fundamental  Exchange  Rules  is  expressed  in  the 
following  result. 

Theorem  1 : B^  is  maintained  as  a nonsingular  matrix  by  the  addition  and 
deletion  of  arcs  if  and  only  if  the  Fundamental  Exchange  Rules  are  applied. 

Proof:  The  master  basis  tree  maintains  a linearly  independent  super- 
structure which,  by  rooting  each  block  (quasi-tree)  of  B^  at  an  ER,  and 
each  slack  arc  of  A in  B,.  at  the  master  root,  assures  that  B, , is  non- 
singular.  Further,  it  is  readily  verified  that  each  operation  prescribed 
by  the  Fundamental  Exchange  Rules  preserves  these  structural  relation- 
ships. The  primary  requirement  is  to  determine  that  none  of  the  possible 
ways  of  modifying  is  overlooked.  We  will  not  trace  the  details  of  a 
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full  itemization  of  cases,  since  they  are  somewhat  tedious,  but  simply 
remark  that  each  of  the  alternatives  can  be  directly  observed  to  be 
handled  correctly  by  the  unifying  construction  of  the  master  basis 
tree.  QED. 

It  may  be  noted  that  Rule  4 of  the  Fundamental  Exchange  Rules,  while 

not  directly  concerned  with  the  addition  and  deletion  of  arcs,  can  affect 

the  density  of  821*  an<*  indirectly  the  composition  of  V.  The  density  of 

any  row  of  B_,  associated  with  A can  be  minimized  (reduced  to  a single 
zl  mn 

non-zero,  if  any  non-zeros  exist)  by  a single  application  of  Rule  4 to 
swap  the  status  of  the  associated  ER  node  with  that  of  the  last  node  of 
its  subtree  (i.e.,  swapping  a row  of  with  a row  of  B^) , using  the 
last  node  function  [5].  Note  that  by  successively  applying  this  pro- 
cedure each  ER  could  be  structured  so  that  it  has  only  one  arc  of  A 

mn 

incident  on  it. 

More  important  and  compelling  is  the  issue  of  whether  there  exists 
an  easily  specifiable  way  to  apply  the  Fundamental  Exchange  Rules  to 
maximize  the  dimension  of  B^  (the  number  of  components  of  xg^)  or  whether, 
due  to  the  influence  of  the  non-network  basis  structure,  a particular 
configuration  for  B^  cannot  be  augmented  except  by  a complex  strategy  of 
removing  current  arcs  and  adding  new  ones.  The  resolution  of  this  issue, 
which  relies  on  the  fact  that  maximizing  the  dimension  of  B^  is  equi- 
valent to  minimizing  the  number  of  EA's  in  the  master  basis  tree,  is  ex- 
pressed in  the  following  result. 

Theorem  2 : The  dimension  of  B^  is  maximized  (and  the  number  of  EA's 
minimized)  by  successively  applying  Rule  1 of  the  Fundamental  Exchange 
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Rules  until  no  more  arcs  are  admissible  to  be  added  by  this  rule. 


Proof : Suppose  Instead  a master  basis  tree  is  obtained  that  cannot  admit 
the  addition  of  any  arc  by  Rule  1,  yet  that  another  master  basis  tree 
exists  with  fewer  EA's.  Since  any  tree  can  be  transformed  into  any 
other  by  some  sequence  of  pairwise  arc  exchanges,  producing  a tree  at 
each  stage,  we  can  find  some  tree  in  such  a sequence  which  allows  no 
swap  reducing  the  number  of  EA's,  and  another  tree  T^,  obtained  by  a 
swap  in  T^,  such  that  a swap  in  will  reduce  the  number  of  EA's.  (Any 
single  swap  that  reduces  the  number  of  EA's  is  an  application  of  Rule  1.) 
Thus,  the  tree  with  a smaller  number  of  EA's  is  two  swaps  away  from  T^. 

But  a basic  result  of  tree  exchanges  is  that  if  two  trees  are  two  swaps 

apart,  then  the  swaps  can  be  executed  in  either  sequence,  producing  a 
tree  at  each  step  (see,  e.g.,  [22]).  This  proves  the  theorem  by  contra- 
diction. QED. 

Theorem  2,  which  applies  either  to  adding  incoming  arcs  to  x^,  or 

transferring  arcs  from  x^  to  also  makes  it  possible  to  maintain 

at  its  maximum  dimension  at  each  Iteration  of  the  primal  simplex  method, 
as  will  be  demonstrated  subsequently.  For  the  special  case  in  which  the 
LP/embedded  network  problem  has  no  non-network  constraints  (though  any 
number  of  non-network  variables),  the  following  observation  is  useful. 


Corollary  1 : If  the  nodes  of  A span  all  the  rows  of  A (i.e.,  if  A^n 
is  empty) , then  every  basic  arc  variable  can  be  Included  in  x^. 

Proof : Any  basis  arc  variable  that  cannot  be  included  in  xR1  must  create 
a cycle  that  does  not  include  an  EA.  This  implies  linear  dependence  in 
the  arc  variable  columns  of  B.  These  columns  must  contain  zeros  in  B^ 
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and  b22*  otherwise  their  arcs  would  intersect  an  ER  and  their  loop  would 
contain  an  EA.  Then  B must  be  singular,  contrary  to  the  fact  that  it 
is  a basis. 

The  value  of  this  Corollary  is  that  it  allows  all  basic  network  arcs 
automatically  to  be  included  in  xfi^  for  LP/embedded  network  problems 
without  side  constraints  (q  = 0),  avoiding  the  work  of  checking  to  see 
whether  the  inclusion  of  any  particular  arc  is  admissible. 

3. 3 Labeling  Algorithms  for  Accelerating  Basis  Computations 

Drawing  on  the  network  topology  of  B,  as  embodied  in  the  construction 
of  the  master  basis  tree,  we  turn  now  to  the  determination  of  special  algo- 
rithms for  processing  this  master  tree.  The  algorithms  are  specifically 
designed  to  carry  out  computations  involving  B ^ that  are  required  in  the 
steps  of  the  primal  simplex  method  (pricing  out  the  basis,  determining 
the  representation  of  the  incoming  variable,  etc.).  In  particular,  we  are 
concerned  with  identifying  the  most  effective  way  to  take  advantage  of 
the  network  structure  of  embedded  in  the  partitioned  inverse. 

The  principal  computations  of  the  primal  simplex  method  that  involve 
B can  be  segregated  into  three  classes  concerned  with  computing: 

(1)  B-^  G,  (2)  HB-^,  and  (3)  HB~*  G.  By  allowing  H and  G in  each  of 
these  classes  to  be  either  a vector  or  a matrix  (and  in  the  case  of  a 
vector,  to  have  either  single  or  multiple  non-zero  components),  it  is 
possible  to  express  the  generic  forms  of  all  of  the  simplex  method  calcu- 
lations involving  B^  in  the  partitioned  inverse.  In  this  section  we 
provide  the  procedures  capable  of  performing  these  calculations  in  a 
manner  that  most  fully  exploits  the  structure  of  B^. 


Significantly,  it 
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turns  out  that  the  most  effective  list  processing  and  labeling  procedures 
differ  in  each  of  the  several  alternative  cases.  The  list  structures  and 
functions  used  in  these  procedures  are  those  commonly  employed  in  network 
optimization.  For  completeness,  we  include  explicit  descriptions  of 
these  lists. 

To  provide  a common  visual  frame  of  reference,  the  root  node  r 
may  be  viewed  as  the  highest  node  in  the  tree  with  all  the  other  nodes 
hanging  below  it.  The  tree  is  then  represented  by  keeping  a pointer  list 
which  contains  for  each  node  w (other  than  the  root)  the  unique  node  v 
above  node  w which  is  connected  to  w by  an  arc.  This  upward  pointer  is 
called  the  predecessor  of  node  w and  will  be  denoted  by  p(w).  Corres- 
pondingly, node  w is  called  an  immediate  successor  of  node  v.  For  con- 
venience, we  will  assume  that  the  predecessor  of  the  root,  p(r),  is  zero. 
Figure  2 illustrates  a tree  rooted  at  node  1,  the  predecessors  of  the 
nodes,  and  other  functions  to  be  described  subsequently.  The  predecessor 
of  a node  is  identified  in  the  p array.  For  example,  the  predecessor  of 
node  16  is  node  5. 

Figure  2 illustrates  additional  tree  information  expressed  in  terms 
of  node  functions,  which  are  often  used  in  computer  implementation  proce- 
dures for  solving  network  problems.  The  first  of  these  functions,  the 
thread  function,  is  denoted  by  t(x).  This  function  is  a "downward"  tree 
pointer.  As  illustrated  in  Figure  2 by  the  dashed  line,  function  t may 
be  thought  of  as  a connecting  link  (thread)  which  passes  through  each 
node  exactly  once  in  a top  to  bottom,  left  to  right  sequence,  starting 
from  the  root  node.  For  example,  in  Figure  2,  t(l)  - 2,  t ( 2 ) - 4, 
t(4)  - 5,  t(5)  - 16,  t (16)  - 8,  etc. 
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Letting  n denote  the  number  of  nodes  In  the  tree,  the  function  t 
satisfies  the  following  inductive  characteristics: 

a)  The  set  {r,t(r),t  (r),...,t  (r)}  is  precisely  the  set  of  nodes 

2 3 2 

of  the  rooted  tree,  where  by  convention  t (r)  - t(t(r)),  t - t(t  (r)), 

k-1 

etc.  The  nodes  r,t(r),...,t  (r)  will  be  called  the  antecedents  of  node 

k,  . 

t (r). 

b)  For  each  node  i other  than  node  tn  ^(r),  t(i)  is  one  of  the 
nodes  such  that  p(t(i))  ■ i,  if  such  nodes  exist.  Otherwise,  let  x de- 
note the  first  node  in  the  predecessor  path  of  i to  the  root  which  has 
an  immediate  successor  y and  y is  not  an  antecedent  of  node  i.  In  this 
case,  t(i)  ■ y. 

c)  tn(r)  * r;  that  is,  the  "last  node"  of  the  tree  threads  back  to 
the  root  node. 

The  reverse  thread  function,  rt(x),  is  simply  a pointer  which  points 
in  the  reverse  order  of  the  thread.  That  is,  if  t(x)  * y,  then  rt(y)  • x. 
Figure  2 also  lists  the  reverse  thread  function  values. 

The  depth  function,  dh(x),  indicates  the  number  of  nodes  in  the  pre- 
decessor path  of  node  x to  the  root,  not  counting  the  root  node  itself. 

If  one  conceives  of  the  nodes  in  the  tree  as  arranged  in  levels  where  the 
root  is  at  level  zero,  and  where  all  nodes  "one  node  away"  from  the  root 
are  at  level  one,  etc. , then  the  depth  function  simply  Indicates  the 
level  of  a node  in  the  tree.  (See  Figure  2.) 

The  cardinality  function,  c(x),  specifies  the  number  of  nodes  con- 
tained in  the  subtree  associated  with  node  x.  By  the  nodes  in  the  sub- 
tree associated  with  node  x,  we  mean  the  set  of  all  nodes  w such  that  the 
predecessor  path  from  w to  the  root  contains  x.  (See  Figure  2.) 
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The  last  subtree  node  function,  f(x),  specifies  the  . „de  in  the 
subtree  of  x that  is  encountered  last  when  traversing  the  nodes  of  this 
subtree  in  "thread  order."  More  precisely,  f(x)  » y where  y is  the 
unique  node  in  the  subtree  of  x such  that  t(y)  is  not  also  a node  in  the 
subtree  of  x.  (See  Figure  2.) 


Note  that  both  the  domain  and  the  range  of  each  of  the  above  dis- 
crete functions  consist  of  a subset  of  the  problem  nodes  and  thus  are 
independent  of  the  number  of  arcs  in  the  LP/embedded  network  problem. 
Since  the  master  basis  tree  contains  m + 1 nodes,  a one-dimensional  array 
of  size  m + 1,  called  a node  length  array , is  allocated  to  each  function 


during  computer  implementation.  The  procedures  for  updating  the  values 
of  the  functions  when  the  tree  is  reconfigured  are  discussed  in  [5,  16]. 


In  computing  a particular  vector  x = B ^G,  the  components  of  x 
are  associated  with  columns  of  via  the  equation  B^x  » G,  and  are 
thus  associated  with  arcs.  Similarly,  the  components  of  the  vector 
W = HB  ^ are  associated  with  rows  of  B^  via  the  equation  WB^  * H and 
are  thus  associated  with  nodes.  Consequently,  in  computing  these  vectors 
we  will  let: 

x^  ■*  the  component  of  x associated  with  the  basis  arc  k (whose  end- 
points are  nodes  k and  p(k)) 

■ the  component  of  W associated  with  node  k. 

These  latter  definitions  will  be  modified  only  to  allow  x and  W to  repre- 
sent matrices  rather  than  vectors , in  which  case  each  row  of  x will  be 
thought  of  as  a row  of  variables  associated  with  the  corresponding  arc, 
and  each  column  of  W as  a column  of  variables  associated  with  the  corres- 
ponding node. 
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Finally,  we  will  refer  to  a basis  arc  as  conformable  if  its  A 

rnn 

direction  agrees  with  its  basis  predecessor  orientation  (i.e.,  the  from 
node  of  the  arc  lies  at  the  predecessor  node  and  the  to  node  of  the  arc 
lies  at  the  successor  node) , and  refer  to  the  arc  as  nonconformable 
otherwise. 


A.  Algorithms  for  computing  x = B nG  (solving  for  x in  the  equation 
B11X  " G^* 

Al.  G is  a column  vector  with  one  non-zero  element.  Let  Gn  equal  the 


non-zero  element  G,  of  G,  associated  with  node  k. 
k 

Al. 1 Let  q ■ p(k).  Let  x^  * GQ  if  the  arc  is  conformable,  and 
let  x^  = -G0  if  the  arc  is  nonconformable. 

Al. 2 If  q is  a root  (either  an  ER  or  the  master  root)  stop.  All 


non-zero  elements  of  x have  been  assigned  the  proper  value.  Other- 


wise, let  k * q and  return  to  Al.l. 


A2.  G is  a column  vector  with  two  non-zero  elements.  Use  the  depth  or 
cardinality  function  in  conjunction  with  the  predecessor  so  that  the 
trace  from  the  two  nodes  that  correspond  to  the  non-zeros  of  G can  be 
interrupted  at  their  intersecting  paths. 


A2.1  Apply  Al  to  each  non-zero,  independently.  Stop  by  the 
criterion  of  A1.2  if  the  paths  do  not  Intersect  first. 

A2.2  If  the  paths  Intersect  before  either  meets  a root,  temporarily 
stop  at  step  Al. 2 when  q is  the  Intersection  node,  and  redefine  G0 
to  be  the  sum  of  the  two  non-zeros  of  G.  If  G0  is  zero,  stop.  Other- 
wise, continue  Al  to  its  termination. 
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A3.  G is  a column  vector  with  more  than  two  non-zero  elements. 

Option  1.  Apply  a generalized  form  of  A2,  not  proceeding  beyond 
any  path  intersection  until  all  paths  meeting  at  that  point  have  been 
traced  to  it.  The  value  of  G0  at  that  point  becomes  the  sum  of  the 
non-zeros  on  the  starting  nodes  of  the  paths  meeting  there.  (This 
can  be  useful  if  most  of  the  paths  are  known  to  Intersect  only  at 
roots. ) 

Option  2.  (Generally  preferred.)  Let  G^  denote  the  i^  element  of 
G (whether  non-zero  or  not).  Using  the  last  node  function,  identify 
the  last  node  of  the  master  tree,  and  designate  this  node  to  be  node 


A3. 1 If  p(k)  is  not  a root,  execute  Step  Al.l  for  GQ  replaced  by 

G.  , and  let  G » G + G . (A3.1  can  be  skipped  if  G,  * 0 thus  calcu- 
k q q k k 

lating  only  non-zero  values.) 

A3. 2 Let  k = rt(k).  If  k is  the  master  root,  stop.  Otherwise, 
return  to  A3.1.  (Note  in  this  procedure,  one  may  avoid  checking 
whether  p(k)  is  a root  in  A3.1  by  allowing  "fictitious"  variables 
x^  to  be  associated  with  EA's.) 

A4.  G is  a matrix  (and  hence  x is  also  a matrix).  Let  G^  denote  the  ith 
row  of  G and  let  x^  denote  the  i^  row  of  x.  Then  apply  algorithm  A3, 
Option  2. 

B.  Algorithms  for  computing  W - HB  (solving  for  W in  the  equation 
WBn  - H.) 


Bl.  H is  a row  vector  with  one  non-zero  element.  Let  H,  denote  the  non- 
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zero  element  of  H,  which  occurs  In  the  kth  position,  associated  with  the 
arc  k of  the  basis  tree  joining  nodes  k and  p(k). 

Bl.l  Let  Hq  * if  the  arc  is  conformable  and  let  Hq  ■ -H^  other- 
wise. Let  i * k. 

B1.2  Let  W - H . 

i o 

Bl. 3 If,  by  the  last  node  function,  node  i is  the  last  node  of  the 
subtree  headed  by  node  k,  stop.  All  non-zeros  of  W have  been  gener- 
ated (with  the  value  H ).  Otherwise,  let  i ■ t(i)  and  return  to 

o 

Bl.  2. 

B2.  H is  a row  vector  with  more  than  one  non-zero  element.  If  the  sub- 
trees containing  the  arcs  associated  with  the  non-zeros  of  H are  known  to 
be  disjoint,  apply  Bl  independently  to  each  subtree.  Otherwise,  select 
any  node  k of  the  master  tree  (possibly  an  ER  or  the  master  root)  which 
heads  a subtree  whose  arcs  "contain"  all  non-zeros  of  H.  Def initionally , 
for  the  following,  let  W ■ 0 if  node  i represents  a root  node  (the  master 
root  or  an  ER) , and  let  H ■ 0 if  arc  i represents  an  externalized  arc 
(which  may  always  be  regarded  as  conformable). 

Option  1.  Let  * 0 and  i ■ t(k).  Let  k*  denote  the  last  node  of 
the  subtree  headed  by  k. 

B2.1  Let  q “ p(i).  Then  let 

W . ■ W + H.  if  arc  i is  conformable 
i q i 

W,  ■ W - H.  if  arc  i is  nonconformable. 
i q i 


2 


B2. 2 Let  i » t(i).  If  i * k*,  stop.  All  non-zero  elements  of  W 
have  been  determined.  Otherwise,  return  to  B2.1. 


Option  2.  (Generally  preferable  if  H has  few  non-zeros.) 

Define  two  lists,  a subtree  header  list  SH(j)  and  a last  node 

list  LN(j),  j = 1,...,J.  To  begin,  let  SH(1)  * k and  let  LN(1)  ■ k*, 

the  last  node  in  the  subtree  headed  by  k.  Let  - 0 and  i ■ t(k). 

B2.3  Let  H = H,  if  arc  i is  conformable  and  let  H ■ -H.  other- 
o i o i 

wise. 

B2.4  Let  W = H . If  i = k*.  go  to  B2.7. 

i o 

B2.5  Let  i = t(i).  If  = 0,  go  to  B2.4.  Otherwise,  let  Hq  “ 

li  + H if  arc  i is  conformable  and  let  H * H - H.  if  arc  i is 
o i o o i 

nonconformable. 

B2.6  If  i * k*,  let  W,  - H and  go  to  B2.7.  Otherwise,  let  SH(J)  ■ i 
i o 

if  the  last  node  in  the  subtree  headed  by  i is  k*  (updating  the  header 
list  to  name  the  most  recent  node  whose  last  node  is  k*) . Otherwise, 
let  J * J + 1,  let  SH(J)  “ i,  let  k*  denote  the  last  node  in  the 
subtree  headed  by  i,  and  go  to  B2.4. 

B2, 7 If  J * 1,  stop.  All  non-zero  components  of  W have  been  deter- 
mined. Otherwise,  let  J * J - 1,  let  k (SH(J),  let  k*  - LN(J)  and 

let  H - W,  . Then  go  to  B2.5. 
ok 

B3.  H is  a matrix  (hence  W is  also  a matrix).  Let  denote  the  i^ 
column  of  H and  let  denote  the  itl1  column  of  W.  Then  apply  algorithm 
B2,  Option  1 or  Option  2.  (If  many  of  the  columns  of  H are  0 vectors. 


1 


Option  2 is  preferable.  If,  further,  non-zero  columns  of  H have  few  non- 


A 
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zero  elements,  then  a variant  of  Option  2 can  be  applied,  particularly 
by  using  pointers  that  name  only  non-zero  elements.) 

Algorithms  Al,  A2,  Bl,  and  B2  Option  1 are  direct  counterparts  of 
algorithms  already  used  in  network  optimization  methods.  Algorithms  A3, 
A4,  B2  Option  2,  and  B3  are  new,  designed  to  handle  the  special  require- 
ments of  the  LP/embedded  network  problem,  with  the  same  types  of  computa- 
tional advantages  that  have  derived  from  their  simpler  predecessors.  One 
of  the  principal  features  shared  by  all  of  these  algorithms,  which  pro- 
vides a primary  basis  for  exploitation  by  the  method  of  the  next  section, 
is  given  in  the  following  result. 

Theorem  3:  The  value  assigned  to  a variable  at  any  iteration  of  any  of 
the  preceding  algorithms  is  the  correct  solution  value  for  the  indicated 
equation  and  is  not  modified  at  any  subsequent  iteration. 

Proof : The  fact  that  a solution  value,  once  assigned,  is  not  changed 
thereafter,  can  be  ascertained  by  tracing  the  steps  of  the  algorithms. 

That  this  value  is  correct,  in  the  case  of  the  new  algorithms,  follows 
from  the  structural  characteristics  of  the  master  basis  tree  already  es- 
tablished, and  from  the  list  processes  employed  in  these  algorithms  (see 
e.g.,  [5]). 

The  usefulness  of  the  "once  and  for  all"  determination  of  values  in- 
dicated in  Theorem  3 manifests  itself  in  the  simplex  SON  optimization  pro- 
cedure by  providing  the  basis  for  three  additional  algorithms  which, 
together  with  the  algorithms  preceding,  yield  the  ability  to  accelerate 
the  calculation  of  more  complex  matrix  products.  The  basis  of  these  algo- 
rithms lies  in  identifying  the  appropriate  manner  to  accumulate  a matrix 
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product  at  Intermediate  stages,  conserving  time  and  memory.  The  form 
of  these  algorithms  is  as  follows. 

C.  Algorithms  for  computing  A ■ HB  G + ZQ 

Cl.  G is  a column  vector  and  H is  a matrix.  Begin  with  Z equal  to  the 


column  vector  ZQ. 


„-l 


Cl. 1 Compute  x - B ^ G by  the  appropriate  member  of  the  A algo- 
rithms (depending  on  the  structure  of  G) . 

Cl. 2 As  each  non-zero  element  x^  of  the  column  vector  X is  computed, 
let  Z * Z 4-  where  is  the  i^  column  of  H. 

C2.  H is  a row  vector  and  G is  a matrix.  Begin  with  Z equal  to  the  row 


vector  ZQ. 


„-l 


C2.1  Compute  W - HB  ^ by  the  appropriate  member  of  the  B algorithms. 
C2.2  As  each  non-zero  element  of  the  row  vector  W is  computed, 
let  A “ A + W^G^,  w^ere  is  the  iC^  row  of  G. 

C3.  G and  H are  both  matrices.  Begin  with  Z equal  to  the  matrix  ZQ. 


Option  1. 


-1 


C3. 1 Compute  x * B ^ G by  algorithm  A4. 

C3. 2 As  each  non-zero  row  x^  of  x is  computed,  let 


Zk  + “ki’i 


Or  alternatively. 


zj  ■ zj  + "ixij 


for  each  row 


of  Z where 


«ki 


is  the 


element  of  the  iC^  column  of  H. 


for  each  column  Z^  of  Z where  x^  is  the 
jth  component  of  x^. 
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Option  2. 

C3.3  Compute  W ■ HB  J by  algorithm  B3. 

C3.4  Aa  each  non-zero  column  of  W is  computed,  let 

+ W , .G.  for  each  row  Z.  of  Z where  W,  . is  the 

ki  i k ki 

kth  element  of  and  is  the  i^ 
row  of  G. 

Or  alternatively, 

Zj  * Zj  + WjGjj  for  each  column  Z^  of  Z where  G_^  is 

the  component  of  G^. 


The  validity  of  these  algorithms  follows  directly  from  the  validity 
of  the  A algorithms  and  B algorithms,  and  from  the  admissible  options  for 
organizing  matrix  computations.  A special  consequence  of  these  methods, 
which  is  uniform  throughout  all  calculations  of  x,  W,  or  Z,  is  as  follows. 

Corollary  2.  The  most  efficient  versions  of  the  A,  B,  and  C algorithms, 
when  either  G or  H is  a matrix,  result  by  storing  G by  row  and  storing 
H by  column. 

Proof:  Immediate  from  the  structure  of  the  algorithms. 

j 

This  outcome  has  noteworthy  implications  for  the  unified  implementa- 
tion of  the  A,  B,  and  C algorithms  in  the  simplex  SON  precedure.  In 
particular,  the  identity  of  G in  the  simplex  SON  procedure  (when  G is  a 

matrix),  is  characteristically  B^»  and  the  Identity  of  H (when  H is  a 

/®2l\ 

matrix),  is  characteristically  B^  or  the  augmented  matrix  \cgi/*  Thus, 
by  Corollary  2,  it  is  preferable  to  store  B^2  by  row  and  B21  by  column. 

This  constitutes  a departure  from  the  methodology  of  previous  compact 

. That  is,  the  result  that  the  most  effective  computa- 


basis  procedures 
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tion  arises  by  storing  different  components  of  B differently,  provides  a 
new  organizational  strategy  for  compact  basis  methods. 

With  these  foundations,  we  are  now  ready  to  specify  the  complete  form 
of  the  simplex  SON  optimization  method. 

4.0  ALGORITHMIC  STEPS 

This  section  describes  in  detail  the  basic  simplex  operations  as 
they  pertain  to  the  LP/embedded  network  problem.  These  operations  in- 
clude initialization,  checking  for  optimality,  finding  the  representation 
of  the  vector  entering  the  basis,  the  basis  exchange  step,  and  calculat- 
ing updated  dual  variable  values.  As  noted  in  Section  3,  it  will  be 
assumed  that  the  basis  is  partitioned  as  in  (11),  that  is  stored  using 
a master  basis  tree,  and  that  V ^ is  the  only  portion  of  the  basis  in- 
verse that  is  being  kept. 

4.1  Initialization 

An  initial  basis  B for  the  LP/embedded  network  problem  can  be  ob- 
tained by  selecting  a set  of  variables  whose  columns  in  A are  linearly 

mn 

independent  and  then  augmenting  these  columns  by  appending  slack  or 
artificial  variables  to  the  problem  that  result  in  satisfying  the  con- 
straints (2)  and  (3).  Columns  for  artificial  variables  whose  unit 
entries  occur  in  the  first  m rows  can  be  treated  as  an  augmentation  of 

the  A matrix, 
mn 

Given  an  initial  basis.  Fundamental  Exchange  Rule  1 can  be  used  to 
partition  the  basis  so  that  the  dimensionality  of  B^  is  maximum. 

Next,  V ■ B22  “ &21  ® 11  ®12  is  comPuted  by  algorithm  C3,  letting 
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Z0  * ®22’  H " -B21*  and  G “ ®12- 


V 1 is,  then,  calculated. 

Once  a basis  has  been  selected  and  partitioned,  the  complementary 
slackness  conditions  can  be  used  to  obtain  dual  variable  values  which 
satisfy: 

B. 


(w  ,w  ) 
m q 


B 


11 

i 

■ 

_L. 

B12 

21 

r 

i 

■ 

■ 

B22  . 

(CB1,CB2) 


(13) 


where  c and  c are  the  vectors  of  objective  function  coefficients  asso- 
BZ 


B1 


XB2* 

In  expanded  form,  (13)  becomes: 

w B 

+ w B»,  ” c_- 

(14) 

m 11 

q 21  Bl 

w B,, 
m 12 

+ Wq®22  “ CB2 

(15) 

Equations  (14)  and  (15)  may  be  rewritten  as  follows: 


WmBll  ’ CB1  * WqB21 

Wq(®22  “ B21B  11B21) 
.-1 


CB2  ~ CB1B  11B12 


Noting  that  v " B22  “ B21B  11B21’  Can  be  stated  as: 


Wq  " (CB2  " CB1B  11B12)V 


(16) 

(17) 

(18) 


w could  be  recomputed  at  future  iterations  using  (18)  and  algo- 

q 

rithm  C2  to  find  -cfi^B  11B12*  Howevert  a8  q (the  number  of  constraints 
in  (3))  increases,  the  number  of  operations  involved  in  this  process 
becomes  prohibitive.  Thus,  for  large  q,  it  is  better  to  treat  w^  as 
an  extra  row  of  V ^ and  update  it  at  each  iteration. 

Given  an  initial  value  of  w^,  (16)  provides  an  excellent  framework 

for  the  calculation  of  w . Once  the  right  hand  side  of  (16)  has  been 

m 

calculated,  w may  be  computed  by  algorithm  B2. 
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4.2  Determination  of  Optimality 


Let  P 


- rp^i 

i Lpq,jJ 


denote  the  j column  vector  of  the  coefficient 


matrix  A,  where  P . is  associated  with  constraints  (2),  and  P . is 

m>J  q.J 

associated  with  constraints  (3).  Optimality  of  both  the  primal  and 
dual  solution  occurs  when  the  dual  solution  (as  well  as  the  primal 
solution)  is  feasible.  A determination  of  this  can  be  made  by  first 
calculating  the  updated  objective  coefficients  c*  associated  with  each 
primal  column  vector  as  follows : 

c*  - c - (w  P + w P ) (19) 

j j ® ®»j  q q»j 


Dual  feasibility  is  achieved  when  the  following  conditions  are  satisfied 


for  all  j : 


cj  > 0,  Xj  - 0 


c*  - 0,  Xj  - basic  (2C 

c»  so,  <J  ■ Y 

If  any  one  of  the  conditions  in  (20)  is  not  satisfied,  then  the  asso- 
ciated primal  column  vector  may  be  selected  to  enter  the  basis. 

4.3  Finding  the  Representation  of  the  Entering  Vector 


Let  P 


8 kJ 


denote  the  column  vector  selected  to  enter  the  basis 


matrix.  The  representation 


a.[“B1' 

aB2. 


of  Pg  in  terms  of  B must  be  com- 


puted so  that  the  vector  to  leave  the  basis  can  be  determined.  This 
computation  involves  solving  the  following  system  of  equations: 

aBl  • »'u  V.  + B'u  “12  V'1  B21  B”n  P„.s  - B'u  “12  V'1  P,,.  (21> 
°B2  ■ V'1<-»21  B'n  P...  + p,.,>-  <2» 


• -r  ‘ ‘v 
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Substituting  into  (21)  yields 

,-l 


aBl  * B 11  (Pm,s  " *12aK2^' 


(23) 


It  is  thus  efficient  to  first  calculate  and  then  to  use  this 
result  to  find  To  compute  first  compute  the  quantity  within 

the  parentheses  of  (22)  by  algorithm  Cl  and  then  multiply  by  V \ To 
compute  a , use  algorithm  A3. 

Using  this  method,  a ratio  test  should  be  performed  immediately  on 
afl2  before  the  computation  of  is  made  since  degenerate  pivots  occur 
frequently  in  network  problems  and,  therefore,  unnecessary  computations 
could  be  avoided. 


4.4  The  Basis  Exchange  Step 

denote  the  vector  selected  to  leave  the  basis.  This 
column  is  identified  by  the  minimum  ratio  calculation.  The  updated  form 
of  this  column,  B ^ Pr,  is  a unit  vector,  whose  non-zero  entry  occurs 
in  the  pivut  row.  To  execute  the  basis  exchange  step,  it  is  necessary 
to  identify  the  segment  of  this  row  in  B * that  affects  the  updating  of 
V \ There  are  two  cases: 

1.  The  outgoing  variable  is  an  element  of  xB2«  The  relevant  seg- 
ment of  the  pivot  row  is  simply  the  row  of  V ^ that  corresponds  to  the 
outgoing  variable. 

2.  The  outgoing  variable  is  an  element  of  xBi.  In  this  case,  the 
pivot  row  segment  that  affects  the  updating  of  V * does  not  lie  in  V ^ 
itself,  but  is  contained  in  the  upper  right  quadrant  of  the  partitioned 
B-*.  Consequently,  from  (12),  this  row  segment  has  the  form 
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-*1  B11  B12  V_1 


(24) 


where  e^  is  the  unit  row  vector  with  a 1 in  the  position  corresponding 


„-l 


to  the  row  of  B associated  with  the  outgoing  variable.  The  computa- 


„-l 


tion  of  (24)  may  be  accomplished  by  first  computing  -e^  B ^ B^»  us*n8 


algorithm  C2  (and  algorithm  B1  in  tie  initial  step).  The  result  is 


then  multiplied  by  V 


-1 


Once  the  appropriate  pivot  row  segment  is  thus  determined,  the 


,.-l 


updating  of  V (and  of  w ) proceeds  in  the  customary  manner  by  a pivot 

q 

step  restricted  to  this  portion  of  B \ 


The  final  type  of  operation  of  the  simplex  SON  method  is  the  modi- 


.-1 


fication  of  V by  the  addition  or  deletion  of  a row  or  column. 


4.5  Changing  the  Dimensionality  of  V 


,-l 


„-l 


The  dimensionality  of  V is  determined  by  the  dimensionality  of  x 


B2 


(hence  also  of  x ),  and  this  in  turn  rests  on  the  application  of  the 
Bl 


Fundamental  Exchange  Rules.  As  previously  noted,  the  initial  composition 
of  and  x^  can  be  determined  by  successive  application  of  Rule  1 to 


maximize  the  dimension  of  B^ , and  consequently  to  minimize  the  dimension 
.-1 


of  V . Once  the  simplex  SON  method  begins  with  this  condition  satisfied. 
Theorems  1 and  2 imply  that  the  basis  exchange  step  can  maintain  this 
condition  at  each  iteration  by  the  transfer  of  at  most  one  element  between 
x^  and  (depending  on  the  configuration  of  the  master  basis  tree, 
which  may  be  modified  by  the  addition  of  an  arc  corresponding  to  the  in- 
coming variable  or  by  the  deletion  of  an  arc  corresponding  to  the  out- 
going variable). 


i 
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In  particular,  by  applying  Rule  1,  2,  or  3 of  the  Fundamental 
Exchange  Rules  to  the  Incoming  and  outgoing  variables  (if  one  or 
both  of  them  correspond  to  arcs) , and  by  applying  Rule  1 to  a 


potential  transfer  variable  (an  element  of  x 0 corresponding  to  an  arc) 

Bz 


— or,  if  appropriate,  applying  Rule  3 to  the  transfer  variable  and  the 
outgoing  variable — the  resulting  reconfiguration  of  the  master  basis 


tree  will  automatically  maximize  the  dimension  of  and  minimize  the 


dimension  of  V \ In  each  case,  V * will  be  modified  by  the  addition 


or  deletion  of  at  most  one  row  and  at  most  one  column  (in  each  combina- 


„-l 


tlon  that  leaves  V square).  The  objective,  then,  is  to  show  how  this 


„-l 


modification  of  V can  be  brought  about. 


„-l 


The  deletion  of  a row  or  column  of  V may  be  accomplished  in  a 
straightforward  manner.  Therefore  we  address  the  operations  involved  in 
the  addition  of  such  a vector.  (In  each  case,  the  additions  should  take 
place  before  the  execution  of  the  pivot  steps,  allowing  the  pivot  step 


to  determine  their  newly  updated  forms.) 


Adding  a row  to  V \ The  only  row  that  may  be  added  to  V ^ is  that 


-1 


of  the  pivot  row,  and  then  only  if  it  does  not  already  lie  in  V (i.e. , 


only  if  the  outgoing  variable  is  an  element  of  x„,).  The  way  to  identify 

B1 


and  calculate  this  updated  row  is  by  (24).  Note  that  the  addition  of  a 


.1 


row  to  V does  not  entail  any  additional  computation,  since  in  this 


case  the  form  of  the  pivot  row  must  be  determined  in  any  event. 


Adding  a column  to  V \ Any  column  of  B ^ not  already  overlapping 


„-l 


V must  lie  in  the  left  half  of  the  partitioned  basis  Inverse  (12). 
Consequently,  the  updated  and  original  form  of  this  column  is  given. 
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r ®bi  i r ek i 

respectively,  by  9 = J and  e^  = |_  q J’ 


where  the  unit  element  of  e. 


and  occur  in  the  same  position.  With  this  identification,  we  may  com- 

pute 6 from  the  explicit  representation  of  the  partitioned  basis  inverse 
(12)  in  precisely  the  manner  used  earlier  to  compute  the  basis  repre- 
sentation a of  the  column  for  the  incoming  variable.  In  particular, 

the  form  of  0 is  obtained  by  substituting  0 , 0 „,  e and  0,  respectively, 

BJL  dZ  k 

for  Qt,..,  0L-,  P , and  P in  (22)  and  (23).  This  yields 

Dl  04  m, s q , s 

6B1  ' B"u  <%  - B12  9B2>  <25> 


6B2  ' V_1  B21  B"u  ek 


The  calculation  of  0^  proceeds  by  first  applying  algorithm  Cl  to  com- 
pute B21  B ^ efc  (and  using  algorithm  A1  in  the  initial  step) , followed 
by  multiplication  by  V \ Then,  unless  both  a row  and  a column  are 
simultaneously  to  be  added  to  V \ it  is  unnecessary  to  compute  0fi^ 
since  0B2  is  the  portion  of  B ^ required. 

On  the  other  hand,  if  a row  is  to  be  added  to  V ^ as  well  as  a 

column,  it  is  necessary  to  compute  the  pivot  row  element  in  0 ..  Thus, 

Bi 

letting  e^  denote  the  unit  row  vector  as  previously  defined  in  the  genera- 
tion of  the  updated  pivot  row,  this  element  of  0 is  given  by 

B1 

el  B"ll  <\  - B12  V <22> 

This  can  be  computed  by  algorithm  C2,  using  algorithm  B1  in  the  first 
part,  and  noting  that  the  matrix  stipulation  for  G may  be  replaced  by  the 
stipulation  that  G is  a column  vector. 


} 
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The  use  of  Che  A,  B,  and  C algorithms  In  these  calculations  materially 
accelerates  the  steps  of  the  staples  SON  method,  following  the  Indicated 
prescriptions. 

5.0  IMPLEMENTATION  AND  CCMPUTA  10NAL  TESTING 

5.1  Iapleaentation 

We  have  implemented  a preliminary  FORTRAN  version  of  the  simpler  SOS 
method  for  capacitated  LP  problems  where  A^  “ 0.  This  tr*-core  code, 
called  PNET/LP,  employs  super-sparsity  (i.e.,  it  stores  only  che  ui'i^u 
non-zero  elements  of  A)  and  keeps  V ^ in  product  form.  'Zii  i*t  also  cr.- 
ploys  the  predecessor,  thread,  reverse  thread,  cartU:»  l tty,  and  lest 
node  functions. 

PNET/LP  has  been  run  on  a CDC  6600,  a CYBER- 74,  and  ar.  « u>AML  V-6. 

The  computer  memory  space  utilized  by  PNET/LP  depends  on  several  factors 

including  (i)  the  computer,  (ii)  the  number  of  unique  non-zeros  in  the 

original  problem,  (iii)  the  number  of  unique  non-zeros  in  the  ETA  file. 

Consequently,  it  is  impossible  to  specify  in  simple  terms  an  exact 

formula  for  the  amount  of  memory  required  by  the  program.  An  approximate 

formula  for  IBM  370  and  AMDAHL  computers  is  46  bytes  per  network  row, 

40  bytes  per  LP  row,  12  bytes  for  each  arc  variable  (i.e.  , for  each 

element  of  x^)  plus  4 bytes  for  each  non-zero  coefficient  in  its  LP 

rows,  and  8 bytes  for  each  non-arc  variable  (i.e.,  for  each  element 

of  x ) plus  4 bytes  for  each  non-zero  coefficient  in  its  LP  rows. 

P 

In  addition,  PNET/LP  keeps  a variably  dimensioned  working  space  array 
which  contains  the  pool  of  unique  non-zero  elements  of  A and  V 1. 


L 


38 


PNET/LP  first  optimally  solves  the  network  portion  of  the  LP / 
embedded  network  problem.  This  optimal  network  basis  is  then  augmented 
by  appropriate  slack  or  artificial  variables  to  form  a starting  basis 
for  the  entire  problem.  During  the  solution  of  the  network  problem, 
PNET/LP  employs  the  modified  row  minimum  start  [19],  and  the  standard 
Phase  I-1I  method  for  handling  artificial  variables.  If  the  optimal 
network  basis  is  augmented  by  artificial  variables,  PNET/LP  minimizes 
the  sum  of  infeasibility  in  Phase  I for  the  full  LP/embedded  network 
problems. 

5.2  Computational  Testing 

In  order  to  evaluate  the  computational  merits  of  PNET/LP,  we  tested 

the  following  three  classes  of  problems:  (i)  pure  network,  (ii)  GUB/LP 

problems,  and  (iii)  embedded  network/LP  problems  where  A ■ 0 and  m is 

mp 

large  relative  to  q. 

The  first  class,  pure  network,  was  selected  in  order  to  determine 
the  relative  efficiency  of  PNET/LP  to  a state-of-the-art  special  purpose 
code  for  solving  network  problems.  To  conduct  our  comparison  we  used 
the  network  code  PNET-I  of  [19]  and  modified  its  pivot  criterion  to 
correspond  to  that  of  PNET/LP.  Our  analysis  disclosed  that  PNET-I  is 
only  twice  as  fast  as  PNET/LP  on  pure  network  problems.  This  is  sur- 
prising in  view  of  the  fact  that  PNET-I  is  150-200  times  faster  than 
commercial  LP  codes  such  as  APEX-III  and  MPSX-370.  However,  since  PNET/LP 
is  designed  to  exploit  network  structures,  it  is  relevant  to  identify 
the  reasons  for  the  difference  in  speed  between  PNET/LP  and  PNET-I. 

These  reasons  are  as  follows: 


. v • - 
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1.  PNET/LP  uses  double  precision  flowing  point  arithmetic  while  PNET-I 
uses  Integer  arithmetic. 

2.  PNET/LP,  due  to  its  use  of  super-sparsity  and  its  ability  to  solve 
any  general  LP  problem,  stores  the  original  problem  data  in  a more  com- 
plex manner  than  PNET-1.  Consequently,  PNET/LP  requires  more  time  to 
access  the  original  problem  during  basis  exchanges  and  pricing  opera- 
tions. 

3.  For  numerical  reasons,  PNET/LP  uses  the  standard  Phase  I-II  method, 
while  PNET-I  uses  the  BIG-M  method.  Computational  testing  [19]  has 
shown  that  the  BIG-M  method  is  substantially  more  efficient  for  pure  net- 
work. problems.  Because  of  these  factors,  the  difference  in  solution  speed 
between  PNET/LP  and  PNET-I  is  much  less  than  might  be  expected.  In 
general,  we  found  that  PNET/LP  is  able  to  solve  network  problems  with  2000 
nodes  and  9000  arcs  in  less  than  15  seconds  using  a FORTRAN  G compiler  on 
an  AMDAHL  V6. 

The  second  class  of  problems,  GUB/LP  problems,  was  selected  because 
generalized  upper  bound  constraints  can  be  viewed  as  a very  simple  form 
of  network  constraints.  Further,  since  the  GUB  feature  has  been  dropped 
from  most  of  the  major  commercial  LP  codes,  we  felt  that  some  evaluation  of 
PNET/LP  on  GUB  problems  would  be  of  interest  to  practitioners.  Our  test 
problems  were  furnished  by  a major  airplane  manufacturer.  The  GUB  portion 
represents  the  assignment  of  plane-types  to  routes.  The  typical  problem 
contained  80  GUB  rows,  14  non-network  rows,  and  130  arc  variables.  On 
these  small  problems  PNET/LP  was  four  times  faster  than  MPSX/370. 

The  third  class  of  problems  is  the  one  which  PNET/LP  is  specifically 
designed  to  solve.  Most  of  our  computational  experience  with  this  problem 
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class  is  for  real  problems  which  were  solved  by  Agrico  Chemical  Company. 
These  problems  involve  the  determination  of  optimal  production  and  dis- 
tribution schedules,  and  were  solved  on  an  AMDAHL  V-6  using  a FORTRAN  G 
compiler.  Unfortunately,  it  was  not  possible  to  compare  PNET/LP  on 
these  problems  to  another  code  because  of  the  proprietary  nature  of 
the  problem  data.  (In  addition,  it  is  difficult  to  obtain  free  computer 
time  on  an  IBM  370/168  or  an  AMDAHL  V-6  for  the  purpose  of  benchmarking 
against  MPSX,  MPSX-370,  or  MPS-III.)  Table  I contains  typical  solution 
statistics  on  Agrico' s three  largest  product  models. 

Subsequent  to  this  testing,  Agrico  acquired  a FORTRAN  H compiler. 
Agrico  has  determined  that  the  FORTRAN  H compiler  reduces  total  run 
time,  including  all  I/O  by  approximately  45%.  Consequently,  the  times  in 
Table  I would  probably  be  reduced  by  45%  using  the  FORTRAN  H compiler. 
Furthermore,  Agrico' s comparison  of  PNET-I  and  PNET/LP  using  the  H compiler 
indicates  that  PNET/LP  is  only  20%  slower  than  PNET-I. 


TABLE  I 

SOLUTION  STATISTICS  ON  PNET/LP 


No.  of  Constraints 

No.  of 

Variables 

PNET/LP* 

No.  of  Network 
Rows  = m 

No.  of  Non-Network 
Rows  = q 

■■ 

Total  Time** 

Total  Pivots 

3179 

20 

15,831 

40 

103  seconds 

10,248 

3442 

6 

21,898 

12 

180  seconds 

17,817 

6192 

10 

21,939 

20 

351  seconds 

10,345 

* AMDAHL  V-6  with  FORTRAN  G compiler 
**  Including  all  I/O  processing 
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In  order  to  provide  some  comparison  of  PNET/LP  to  a commercial  LP  code 
on  the  third  class  of  problems,  we  solved  one  randomly  generated  problem  on 
the  CYBER-74  computer  using  PNET/LP  and  APEX-III.  Table  II  contains  the 
problem  specifications  and  indicates  that  this  problem  can  be  solved  at  least 
70  times  less  expensively  with  PNET/LP  than  with  APEX-III. 

TABLE  II 

PNET/LP  VS.  APEX-III 


No.  of 

Constraints 

No.  of 

Variables 

PNET/LP 

APEX-III 

No.  of  Network 
Rows  = m 

No.  of  Non-Network 
Rows  = q 

tm 

1000 

1 

5000 

i 

$3.11* 

$210.68** 

* AMDAHL  V-6  with  FORTRAN  G compiler 
**  Including  all  I/O  processing 

While  the  above  computational  testing  is  quite  limited,  it  indicates 
that  the  simplex  SON  algorithm  may  be  extremely  efficient  for  solving  large 
embedded  network/LP  problems.  At  this  point,  however,  we  would  caution  the 
reader  that  it  would  be  premature  to  extrapolate  these  results  to  other 
problems  and  in  particular,  to  other  problem  classes.  An  exhaustive  com- 
putational study  is  required  before  any  general  Inferences  can  be  drawn. 
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